function [energy,sumC,sumJ,sumV]=hrcorrcp3(element_name,e2,e1,phi)

% hrcorrcp: return Holm-Ribberfors corrected HF Compton profiles as
% tabulated by Biggs et al.
%INPUT: element_name = name of the element
%                 e2 = analyzer energy [keV]
%                 e1 = incident energy [keV]
%                phi = scattering angle [deg two theta]
%OUTPUT: energy = energy
%        sumC   = HF core compton profile
%        sumJ   = HF compton profile
%        sumV   = HF valence profile
%depends on PZ data in folder 'directory'

directory ='../helpers/pzdata/';

switch upper(element_name)
    case 'C'
        shells    = [1 2 3];
        electrons = [2 2 2];
        valence   = [0 2 2];     
        edges     = [284.2 18 7.2];
        element_name = 'C';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);  
    case 'N'
        shells    = [1 2 3 3];
        electrons = [2 2 2 1];
        valence   = [0 0 2 1];     
        edges     = [409.9 37.3 17.5 17.5];
        element_name = 'N';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);  
    case 'BA'
        shells    = [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18];
        electrons = [ 2 2 2 4 2 2 4 4 6 2 2 4 4 6 2 2 4 2 ];
        valence   = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 4 2 ];
        edges     = [ 37441 5989 5624 5247 1293 1137 1063 795.7 780.5 253.5 192 178.6 92.6 89.9 30.3 17 14.8  0 ];
        element_name   = 'Ba';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);    
    case 'SI'
        shells    = [ 1 2 3 3 4 5 ];
        electrons = [ 2 2 2 4 2 2 ];
        valence   = [ 0 0 0 0 2 2 ];  
        edges     = [ 1839 149.7 99.2 99.8 8 3 ];
        element_name   = 'SI';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);
    case 'O'
        shells    = [ 1 2 3 3 ];
        electrons = [ 2 2 2 2 ]; 
        valence   = [ 0 0 2 2 ]; 
        edges     = [ 543 41.6 18.2 18.7 ];
        element_name   = 'O';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);
    case 'BE'
        shells    = [1 2];
        electrons = [2 2 ];
        valence   = [0 2];
        edges     = [111.5 8];
        element_name = 'BE';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);
    case 'H'
        shells    = [1 ];
        electrons = [1 ];
        valence   = [1 ];
        edges     = [13.6];
        element_name = 'H';
        theory=load([directory 'PZ',upper(element_name),'.DAT']);
    otherwise
        disp('WARNING: element not yet in library!')
        return
end
%momentum transfers and pz axis
q=momtrans2(e1,e2,phi);q=q(:);
pz=e2pz(e1,e2,phi);pz=pz(:);

%interpolate and mirror HF data on pz scale, multiply by no. of electrons
pzhf = theory(:,1);
jhf_pre  = theory(:,1+shells);
jhf = zeros(length(pz),length(shells));
for n=1:length(shells)
    jhf_pre(:,n) = jhf_pre(:,n).*electrons(n);
    jhf(:,n)=interp1(pzhf(:,1),jhf_pre(:,n),abs(pz));
end

% hydrogenic CP after Holm and Ribberfors for completely occupied shells
%gamma 1s
fitfct = @(a)(abs(max(jhf(:,1))-max(8*a^5/3/pi./(a^2+pz.^2).^3)));
gamma1s = fminsearch(fitfct,sum(electrons));
j0(:,1) = 8*gamma1s^5/3/pi./(gamma1s^2+pz.^2).^3;
j1(:,1) = 2*gamma1s*atan2(pz,gamma1s)-3/2*pz; 
j1(:,1)=j1(:,1)./q.*j0(:,1);
%gamma 2s
if length(electrons(1))>1
    fitfct = @(a)(abs(max(jhf(:,2))-max((a^4 - 10*a^2*pz.^2 + 40*pz.^4).*128*a^5/15/pi./(a^2 + 4*pz.^2).^5)));
    gamma2s = fminsearch(fitfct,sum(electrons)*2/3);
    j0(:,2) = (gamma2s^4 - 10*gamma2s^2*pz.^2 + 40*pz.^4).*128*gamma2s^5/15/pi./(gamma2s^2 + 4*pz.^2).^5;
    j1(:,2) = 2*gamma2s*atan2(2*pz,gamma2s)-5/4*(gamma2s^4+48*pz.^4)./(gamma2s^4-10*gamma2s^2*pz.^2+40*pz.^4).*pz; 
    j1(:,2) = j1(:,2)./q.*j0(:,2);
end
%gamma 2p
if length(electrons)>2
    fitfct = @(a)(abs(max(sum(jhf(:,shells(:)==3),2))-max((a^2+20*pz.^2)*64*a^7/5/pi./(a^2+4*pz.^2).^5)));
    gamma2p = fminsearch(fitfct,sum(electrons)*1/3);
    j0(:,3) = (gamma2p^2+20*pz.^2)*64*gamma2p^7/5/pi./(gamma2p^2+4*pz.^2).^5;
    j1(:,3) = 2*gamma2p*atan2(2*pz,gamma2p)-2/3.*pz.*(10*gamma2p^2+60*pz.^2)./(gamma2p^2+20*pz.^2);
    j1(:,3)=j1(:,3)./q.*j0(:,3);

end

%consider not-filled 2p shells
if length(shells)>3
    j0(:,3) = j0(:,3)*electrons(3)/(electrons(3)+electrons(4));
    j0(:,4) = j0(:,3)*electrons(4)/(electrons(3)+electrons(4));
    j1(:,3) = j1(:,3)*electrons(3)/(electrons(3)+electrons(4));
    j1(:,4) = j1(:,3)*electrons(4)/(electrons(3)+electrons(4));
end

%convert to enery scale
energy = (e1-e2)*1000;

%assign profiles, add up j0 and j1 (electron count in 2p shell already accounted for)
HR = j0 +j1;
if length(shells)<5
    J = jhf;
else
    J = [HR jhf(:,5:end)];
end

% discard data below zero
J(energy<0,:)=[];
energy(energy<0)=[];

% cut at the edge
for n=1:length(shells)
  dmy = J(:,n);
  dmy(energy<=edges(n))=0;
  J(:,n) = dmy;
end

% add profiles
sumJ = zeros(length(energy),1);
for n=1:length(shells)
  sumJ = sumJ + J(:,n);
end
% add only valence
sumV = zeros(length(energy),1);
for n=1:length(shells)
  sumV = sumV + valence(n)*J(:,n)/electrons(n);
end
% add only core
sumC = sumJ - sumV;

